#ifndef xpack_math_sin
#define xpack_math_sin
#pragma push_macro("xuser")
#undef  xuser
#define xuser mixc::math_sin::inc
#include"define/base_type.hpp"
#include"math/abs.hpp"
#include"macro/xalign.hpp"
#include"macro/xexport.hpp"
#include"math/literal.hpp"
#include"math/private/parallel_lut.hpp"
#include"utils/array.hpp"
#pragma pop_macro("xuser")

namespace mixc::math_sin{
    template<class type_t, uxx multiple>
    inline xalign(32) auto lut = inc::parallel_lut<multiple>(
        type_t(+1.00000000000000000e+00),
        type_t(-1.66666666666666657e-01),
        type_t(+8.33333333333333322e-03),
        type_t(-1.98412698412698413e-04),
        type_t(+2.75573192239858925e-06),
        type_t(-2.50521083854417202e-08),
        type_t(+1.60590438368216133e-10),
        type_t(-7.64716373181981641e-13),
        type_t(+2.81145725434552060e-15),
        type_t(-8.22063524662432950e-18)
    );

    namespace inc{
        using namespace literal;
    }
}

namespace mixc::math_sin::origin{
    // 注意：
    // 此 unsafe 函数不带定义域检查，需要使用者保证 x 的范围属于 [-pi/2, +pi/2] 的闭区间
    template<class type_t, class ... args_t>
    inline auto sin_unsafe(type_t const & x0, args_t const & ... xn) {
        constexpr auto multiple = 1 + sizeof...(args_t);
        constexpr auto step     = sizeof(type_t) == sizeof(f32) ? 6 : 10;

        auto && x2      = inc::make_array((type_t &)x0, (args_t &)xn...);
        auto && xm      = inc::make_array((type_t &)x0, (args_t &)xn...);
        auto && rx      = inc::array<type_t, multiple>();

        for(uxx i = 0; i < multiple; i++){
            x2[i]      *= x2[i];
        }

        for(uxx i = 0; i < step; i++){
            for(uxx j = 0; j < multiple; j++){
                rx[j]  += xm[j] * lut<type_t, multiple>[i * multiple + j];
                xm[j]  *= x2[j];
            }
        }

        if constexpr (multiple == 1){
            return rx[0];
        }
        else{
            return rx;
        }
    }

    template<class type_t, class ... args_t>
    inline auto sin(type_t const & x0, args_t const & ... xn) {
        auto normalize = [](type_t x){
            constexpr auto hpi  = type_t(inc::pi * 0.5); // half pi
            constexpr auto rhpi = type_t(1.0) / (hpi);

            auto is_neg         = x < 0;

            // 先
            if (inc::abs(x) <= hpi){
                return x; // sin 不关于原点对称，所以不能返回 abs(x)
            }

            // 后
            if (is_neg){
                x               = -x;
            }

            auto q              = i64(x * rhpi);         // 整数倍
            auto r              = type_t(x - q * hpi);   // 余数

            // sin 函数是一个周期函数
            // 一共可以分成 4 个段：
            // 1. [-pi  , -pi/2)
            // 2. [-pi/2,     0)
            // 3. [0,     +pi/2)
            // 4. [+pi/2, +pi  )
            // 其中：
            // 段 1,2 和 段 3,4是关于原点对称的，所以有 -sin(x) = sin(-x) (x > 0)
            // 所以有矫正 2 中判断当前 x 是否在 1,2段
            // 此外：
            // 段 1 和段 2 关于 x = -pi/2 对称
            // 段 3 和段 4 关于 x = +pi/2 对称
            // 所以有 sin(pi/2 - x) = sin(pi/2 + x)  (x >= 0 and x <= pi/2)

            // 矫正1：
            r                   = q & 1 ? hpi - r : r;

            // if (is_neg){
            //     if ((q & 0x2) == 0){
            //         r           = -r;
            //     }
            // }
            // else {
            //     if (q & 0x2){
            //         r           = -r;
            //     }
            // }
            // 
            // 矫正2：
            // 等效上述代码
            if (is_neg ^ ((q & 0x2) != 0)){
                r           = -r;
            }
            return r;
        };
        return sin_unsafe(normalize(x0), normalize(xn)...);
    }

    inline f32 sin_unsafe(f32 x) {
        return sin_unsafe<f32>(x);
    }

    inline f64 sin_unsafe(f64 x) {
        return sin_unsafe<f64>(x);
    }

    inline f32 sin(f32 x) {
        return sin<f32>(x);
    }

    inline f64 sin(f64 x) {
        return sin<f64>(x);
    }
}

#endif

xexport_space(mixc::math_sin::origin)